body .main-container {
  max-width: 100% !important;
  width: 100% !important;
}
body {
  max-width: 100% !important;
}

body, td {
   font-size: 16px;
}
code.r{
  font-size: 14px;
}
pre {
  font-size: 14px
}
knitr::opts_chunk$set(message = FALSE, 
                      warning = FALSE,
                      fig.align = "center",
                      out.width = "100%"
                      )
set.seed(1982)
library(INLA)
# Load the data
data(Seeds)

# Create a data frame
df <- data.frame(y = Seeds$r, Ntrials = Seeds$n, Seeds[, 3:5])

# explore the data
df |> head() |> knitr::kable(caption = "Data")
Data
y Ntrials x1 x2 plate
10 39 0 0 1
23 62 0 0 2
23 81 0 0 3
26 51 0 0 4
17 39 0 0 5
5 6 0 1 6

0.1 Model specification

The model specification is

\[ y_i\sim\text{Binomial}\left(N_i, p_i\right)\label{eq1}\tag{A} \]

\[ p_i = \text{invlogit}(\eta_i) = \frac{\exp(\eta_i)}{1+\exp(\eta_i)} = \frac{1}{1+\exp(-\eta_i)}\qquad \left(\text{i.e., }\eta_i = \text{logit}(p_i) =\log\left(\frac{p_i}{1-p_i}\right) \right) \]

\[ \eta_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} +v_i \]

\[ \beta_i \sim \text{N}(0, 1000) \]

\[ v_i\sim\text{N}(0,\sigma^2) \] \[ \sigma\sim\text{Exp}(\lambda) \qquad \left(\text{i.e., }\pi(\sigma) = \lambda e^{-\lambda\sigma}\right) \]

\[ \lambda \text{ is such that } \pi(\sigma>1) = 0.01\qquad\left(\text{i.e., }0.01 = \pi(\sigma>1) = 1- \pi(\sigma\leq1) = e^{-\lambda}\implies \lambda = -\dfrac{\log(0.01)}{1}\approx4.6\right)\label{eq2}\tag{B} \]

Expressions from \(\eqref{eq1}\) to \(\eqref{eq2}\) define the model.

  • \(v_i\) are the random effects.
  • \(\pi(\sigma>1) = 0.01\) corresponds exactly to list(theta = list(prior = "pc.prec", param = c(1,0.01)))
  • \(x_{1i}\) corresponds to x1
  • \(x_{2i}\) corresponds to x2
  • \(v_i\) corresponds to f()
  • model = "iid" tells inla() we are using random effects
  • INLA does not deal with \(\sigma^2\) but with \(\tau = \sigma^{-2}\). Therefore, using the change of variable formula, we can see that \(\pi(\tau) = \dfrac{\lambda}{2}\tau^{-3/2}e^{-\lambda \tau^{-1/2}}\). Of course, \(\pi(\sigma) = \pi(\tau)\left|\dfrac{d\tau}{d\sigma}\right| = \lambda e^{-\lambda\sigma}\) and we define \(x=\log(\tau)\), then \(\pi(x) = \pi(\tau)\left|\dfrac{d\tau}{dx}\right| = \dfrac{\lambda}{2}\exp\left(-\dfrac{x}{2} - \lambda \exp\left(-\dfrac{x}{2}\right)\right)\)
  • variable plate contains the identity of each observations. That is, plate = 1:nrow(df).

The logit() function is not referenced as "logit" in control.family = list(control.link = list(model = "logit")). This "logit" is an option from INLA.

logit <- function(p){log(p/(1-p))}
invlogit <- function(x){exp(x)/(1 + exp(x))}
# call to inla
formula1 = y ~ x1 + x2 + f(plate, 
                           model = "iid",  
                           hyper = list(theta = list(prior = "pc.prec",
                                                     param = c(1,0.01))))
res1 = inla(formula = formula1, 
            data = df, 
            family = "binomial",
            Ntrials = Ntrials, 
            control.family = list(control.link = list(model = "logit")),
            control.predictor = list(compute = TRUE, link = 1))

0.2 SUMMARY

Here we show how to read the summary

#str(res1, 1)
summary(res1)
## Time used:
##     Pre = 0.375, Running = 0.155, Post = 0.0129, Total = 0.543 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.387 0.181     -0.739   -0.390     -0.017 -0.390   0
## x1          -0.360 0.231     -0.838   -0.353      0.076 -0.353   0
## x2           1.033 0.221      0.589    1.034      1.469  1.034   0
## 
## Random effects:
##   Name     Model
##     plate IID model
## 
## Model hyperparameters:
##                      mean    sd 0.025quant 0.5quant 0.975quant mode
## Precision for plate 19.78 32.61       2.96    10.55      84.94 6.06
## 
## Marginal log-Likelihood:  -68.44 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
res1$summary.random$plate |> head() |> knitr::kable(format = "markdown", caption = "Summary for the random effects $v_i$. Keep in mind that each of the $v_i$ has a marginal posterior distribution")
Summary for the random effects \(v_i\). Keep in mind that each of the \(v_i\) has a marginal posterior distribution
ID mean sd 0.025quant 0.5quant 0.975quant mode kld
1 -0.3048725 0.2799168 -0.9193189 -0.2777256 0.1690210 -0.1865579 1.8e-06
2 -0.0859021 0.2322271 -0.5734249 -0.0753286 0.3594508 -0.0374485 2.0e-07
3 -0.3302246 0.2477875 -0.8626290 -0.3116615 0.0969325 -0.2731563 1.2e-06
4 0.2247962 0.2451080 -0.2194151 0.2075764 0.7478825 0.1644008 6.0e-07
5 0.0559880 0.2459764 -0.4318958 0.0503295 0.5620322 0.0481711 2.0e-07
6 0.1082030 0.3233224 -0.5018691 0.0848421 0.8176969 0.0357540 1.6e-06 
res1$summary.fitted.values |> head() |> knitr::kable(format = "markdown", caption = "Summary for the fitted values $p_i$. Keep in mind that each of the $p_i$ has a marginal posterior distribution")
Summary for the fitted values \(p_i\). Keep in mind that each of the \(p_i\) has a marginal posterior distribution
mean sd 0.025quant 0.5quant 0.975quant mode
fitted.Predictor.01 0.3366888 0.0585931 0.2213399 0.3376109 0.4489894 0.3417845
fitted.Predictor.02 0.3851974 0.0491712 0.2903125 0.3845965 0.4840581 0.3838758
fitted.Predictor.03 0.3300258 0.0470660 0.2394589 0.3297444 0.4221690 0.3291541
fitted.Predictor.04 0.4597625 0.0579002 0.3544738 0.4570430 0.5794844 0.4501589
fitted.Predictor.05 0.4190118 0.0580471 0.3108036 0.4165443 0.5398323 0.4108377
fitted.Predictor.06 0.6748795 0.0746330 0.5208798 0.6752715 0.8202451 0.6717140 
res1$summary.linear.predictor |> head() |> knitr::kable(format = "markdown", caption = "Summary for the linear predictors $\\eta_i$. Keep in mind that each of the $\\eta_i$ has a marginal posterior distribution")
Summary for the linear predictors \(\eta_i\). Keep in mind that each of the \(\eta_i\) has a marginal posterior distribution
mean sd 0.025quant 0.5quant 0.975quant mode kld
Predictor.01 -0.6900710 0.2694953 -1.2578753 -0.6739590 -0.2047550 -0.6301495 7e-07
Predictor.02 -0.4725910 0.2100300 -0.8938666 -0.4700832 -0.0637893 -0.4697858 1e-07
Predictor.03 -0.7160160 0.2160388 -1.1556482 -0.7093413 -0.3138757 -0.7088185 3e-07
Predictor.04 -0.1633026 0.2358492 -0.5994313 -0.1722528 0.3206572 -0.2093541 3e-07
Predictor.05 -0.3313657 0.2414093 -0.7963651 -0.3369757 0.1596674 -0.3380490 2e-07
Predictor.06 0.7522351 0.3550642 0.0835679 0.7321254 1.5180090 0.6785054 1e-07 
data.frame(p = res1$summary.fitted.values$mean, 
           eta = res1$summary.linear.predictor$mean, 
           p.from.eta = invlogit(res1$summary.linear.predictor$mean),
           eta.from.p = logit(res1$summary.fitted.values$mean)) |>
  head() |>
  knitr::kable(format = "markdown", caption = "This shows that `res1$summary.fitted.values` corresponds to the variable ($p_i$) related to the linear predictor ($\\eta_i$) via the link function (`logit()`), and `res1$summary.linear.predictor` corresponds to the linear predictor ($\\eta_i$)")
This shows that res1$summary.fitted.values corresponds to the variable (\(p_i\)) related to the linear predictor (\(\eta_i\)) via the link function (logit()), and res1$summary.linear.predictor corresponds to the linear predictor (\(\eta_i\))
p eta p.from.eta eta.from.p
0.3366888 -0.6900710 0.3340173 -0.6780850
0.3851974 -0.4725910 0.3840032 -0.4675454
0.3300258 -0.7160160 0.3282709 -0.7080684
0.4597625 -0.1633026 0.4592648 -0.1612990
0.4190118 -0.3313657 0.4179084 -0.3268315
0.6748795 0.7522351 0.6796655 0.7303383 
# betas
res1$summary.fixed |> head() |> knitr::kable(format = "markdown", caption = "Summary for the fixed effects")
Summary for the fixed effects
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.3869998 0.1812751 -0.7387982 -0.3901860 -0.0174492 -0.3900842 0e+00
x1 -0.3598794 0.2306549 -0.8379333 -0.3525628 0.0760315 -0.3529318 1e-07
x2 1.0328605 0.2214283 0.5891781 1.0343050 1.4690759 1.0344034 0e+00 
# hyperparameters
rbind(res1$summary.hyperpar, res1$internal.summary.hyperpar, "exp(Log precision for plate)" = exp(res1$internal.summary.hyperpar)) |> 
  knitr::kable(format = "markdown", caption = "Summary for the hyperparameters")
Summary for the hyperparameters
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for plate 19.776422 32.6088461 2.962334 10.552322 84.934835 6.061489
Log precision for plate 2.453651 0.8231339 1.085977 2.356346 4.441884 2.167514
exp(Log precision for plate) 11.630734 2.2776265 2.962334 10.552322 84.934835 8.736543

This is how the output of inla() corresponds to our model specification.

  • \(p_i\) = res1$summary.fitted.values (just the summary, not the whole the marginal posterior)

  • \(\eta_i\) = res1$summary.linear.predictor (just the summary, not the whole the marginal posterior)

  • \(\beta_i\) = res1$summary.fixed (just the summary, not the whole the marginal posterior)

  • \(\tau\) = res1$summary.hyperpar (just the summary, not the whole the marginal posterior)

0.3 MARGINAL POSTERIOR

0.3.1 For the fixed effects

Let us get the marginal posterior for \(\beta_1\). No transformation is needed in this case.

marginal.posterior.beta1.few = res1$marginals.fixed$x1 # we just get a few
marginal.posterior.beta1.all = inla.tmarginal(fun = function(x) x, marginal = res1$marginals.fixed$x1) # we get more
# let us check
data.frame(mean.from.summary = res1$summary.fixed$mean[2], 
           mean.from.marginal.posterior.few = mean(marginal.posterior.beta1.few[,1]),
           mean.from.marginal.posterior.all = mean(marginal.posterior.beta1.all[,1])) |> t() |>
  knitr::kable(format = "markdown", caption = "Checking that we are getting the marginal posterior correctly")
Checking that we are getting the marginal posterior correctly
mean.from.summary -0.3598794
mean.from.marginal.posterior.few -0.3670172
mean.from.marginal.posterior.all -0.3602327 
plot(marginal.posterior.beta1.all, type = "p", xlab = expression(beta[1]), ylab = "Probability density", col = "red")
points(marginal.posterior.beta1.few, col = "blue")
legend("topleft", pch = 1, col = c("red", "blue"), legend = c("all samples", "few samples"))

To get the marginal posterior in a meaningful scale, we should change fun from inla.tmarginal(fun, marginal) as appropriate.

0.3.2 For the hyperparameters

Recall that inla() works with precision and that the internal parameters are in log() scale. So if we want the marginal posterior for \(\sigma\) we need to use a transformation as appropriate.

Recall that \(\tau = \sigma^{-2}\).

So if we let \(\tau\) = res1$marginals.hyperpar$Precision for plate, then to get \(\sigma = \tau^{-1/2}\), we need to use the transformation fun(x) = x^(-1/2).

If we let \(\log(\tau)\) = res1$internal.marginals.hyperpar$Log precision for plate, then to get \(\sigma = \tau^{-1/2}\), we need to use the transformation fun(x) = e^(-x/2).

# posterior of sigma
marginal.posterior.sigma.from.tau = inla.tmarginal(fun = function(x) x^(-0.5), marginal = res1$marginals.hyperpar$`Precision for plate`)
marginal.posterior.sigma.from.log.tau = inla.tmarginal(fun = function(x) exp(-x/2), marginal = res1$internal.marginals.hyperpar$`Log precision for plate`)

plot(marginal.posterior.sigma.from.tau, type = "p", xlab = expression(sigma[iid]), ylab = "Probability density", col = "red")
points(marginal.posterior.sigma.from.log.tau, col = "blue")
legend("topright", pch = 1, col = c("red", "blue"), legend = c("sigma from tau", "sigma from log tau"))

So if we let \(\tau\) = res1$marginals.hyperpar$Precision for plate, then to get \(\tau\), we need to use the transformation fun(x) = x.

If we let \(\log(\tau)\) = res1$internal.marginals.hyperpar$Log precision for plate, then to get \(\tau\), we need to use the transformation fun(x) = e^x.

# posterior of tau
marginal.posterior.tau.from.prec = inla.tmarginal(fun = function(x) x, marginal = res1$marginals.hyperpar$`Precision for plate`)
marginal.posterior.tau.from.log.prec = inla.tmarginal(fun = function(x) exp(x), marginal = res1$internal.marginals.hyperpar$`Log precision for plate`)

plot(marginal.posterior.tau.from.prec, type = "p", xlab = expression(tau), ylab = "Probability density", col = "red")
points(marginal.posterior.tau.from.log.prec, col = "blue")
legend("topright", pch = 1, col = c("red", "blue"), legend = c("tau from prec", "tau from log prec"))

mean(marginal.posterior.tau.from.prec[,1])
## [1] 17.96846
mean(marginal.posterior.tau.from.log.prec[,1])
## [1] 17.99411
mean(res1$marginals.hyperpar$`Precision for plate`[,1])
## [1] 27.84323
mean(exp(res1$internal.marginals.hyperpar$`Log precision for plate`[,1]))
## [1] 27.84323

0.3.3 For the random effects

m.plate.7 = inla.tmarginal(fun = function(x) x, marginal = res1$marginals.random$plate$index.7)
plot(m.plate.7, type="l", xlab = "marginal plate nr 7", ylab = "Probability density")

Here are different ways to get \(\eta_7 = \beta_0 + \beta_1x_{17} + \beta_2x_{27} + v_7\)

res1$summary.linear.predictor$mean[7]
## [1] 0.8104082
res1$summary.fixed$mean[1] + res1$summary.fixed$mean[2]*df[7,3] + res1$summary.fixed$mean[3]*df[7,4] + mean(res1$marginals.random$plate$index.7[,1])
## [1] 0.8140505
res1$summary.fixed$mean[1] + res1$summary.fixed$mean[2]*df[7,3] + res1$summary.fixed$mean[3]*df[7,4] + res1$summary.random$plate$mean[7]
## [1] 0.8111264